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Abstract: Omitted variable bias can affect treatment effect estimates ob- 
tained from observational data due to the lack of random assignment to 
treatment groups. Sensitivity analyses adjust these estimates to quantify 
the impact of potential omitted variables. This paper presents methods 
of sensitivity analysis to adjust interval estimates of treatment effect — 
both the point estimate and standard error — obtained using multiple lin- 
ear regression. Central to our approach is what we term benchmarking, 
the use of data to establish reference points for speculation about omit- 
ted confounders. The method adapts to treatment effects that may differ 
by subgroup, to scenarios involving omission of multiple variables, and to 
combinations of covariance adjustment with propensity score stratification. 
We illustrate it using data from an influential study of health outcomes of 
patients admitted to critical care. 
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1. Introduction 

1.1. Methodological context 

In a common use of multiple linear regression, one regresses an outcome variable 
on an treatment variable and adjustment variables, then interprets the fitted 
treatment-variable coefficient as an estimate of the treatment's effect on the 
outcome. The interpretation relies on the assumptions of the linear model and 
some assumption to the effect that there either are no unmeasured confounders 
or at least none that demand adjustment. (Ignorability of treatment assignment 
(Rubin, 1978; Holland, 1988, Appendix) is one such assumption; there are many 
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variants.) The linearity assumptions are often testable given the data, but the 
remaining assumption is not. When regression results are questioned, it's often 
the non-confounding assumption that is the focus of doubt. 

Because the issue arises even with the most thorough observational studies, 
adjusting for any number of covariates, it fuels cynicism about observational re- 
search. If the possibility of unmeasured variable bias can't be removed, then why 
bother with potential confoundcrs, particularly those that are difficult to mea- 
sure, or not obvious threats? It might be clear that the damage from omitting 
a confounder W would be reduced by adjustment for available correlates of W 
yet, because introducing these correlates would draw attention to the absence of 
W, not at all clear that effecting the additional adjustments would enhance the 
credibility of the research. Plainly, the problem here is not the methodological 
strategy of broadly adjusting for relevant baseline characteristics but an absence 
of, or lack of awareness of, suitable methods with which to quantify benefits of 
more comprehensive confounder controls. 

Sensitivity analyses, procedures in some way quantifying the degree of omit- 
ted variable bias needed to nullify or reverse key conclusions of a study, should be 
adaptable to this purpose. Sensitivity analysis methods for various models and 
data structures are proposed in Cornfield et al. (1959), Rosenbaum and Rubin 
(1983), Rosenbaum (1988), Copas and Li (1997), Robins et al. (2000), Scharfstcin and Irizarry 
(2003), Marcus (1997), Lin et al. (1998), Frank (2000) andlmbens (2003), among 
others, the last four bearing closest resemblance to the approach to be presented 
here. Invariably the methods start by in some way quantifying relationships 
between hypothetical omitted variables and included ones, go on to give an al- 
gorithm for converting these parameters into impacts on estimates, p- values or 
confidence limits, and then leave to researchers themselves the task of decid- 
ing what parameter values are plausible or relevant. Here we develop a method 
following the first parts of this general recipe but then doing a bit more to 
help researchers calibrate intuitions about speculation parameters. The result- 
ing technique possesses a combination of advantages making it both uniquely 
practicable and conducive to certain insights. First, it applies to inferences made 
with ordinary multiple regression, as we show in Section 2, as well as to infer- 
ences made with regression in combination with propensity-score stratification, 
a topic discussed in Section 4.3. Second, it quantifies relationships between omit- 
ted and included variables in terms intrinsic to multiple regression, permitting 
intuitions for the relationships to be calibrated with a few additional regression 
fits (§ 2). (Angrist and Imbens (2002, §4) suggest such calibration in a related 
context.) Third, they represent effects of the omission with just two quantities, 
one tracking confoundedness with the treatment variable and the other measur- 
ing conditional association with the outcome (§ 2). Fourth, there are important 
practical benefits to sensitivity analysis based on both of these quantities - 
"dual" sensitivity analyses, in Gastwirth et al.'s terminology (1998) — in that 
analysis based only on confoundedness with the treatment variable is likely to 
overstate sensitivity. Our method makes this plain (§ 3), as we will demonstrate 
with a case study to be introduced presently, although some other methods may 
obscure it. Fifth, they give closed-form characterizations of how confidence in- 
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tervals, as opposed only to estimates or hypothesis tests, could be changed by 
inclusion of the omitted confounder (§ 3.2). Sixth, the method readily adapts 
to analyses in which several omissions are suspected, or where interactions with 
the treatment variable are used to handle possible effect heterogeneity (§ 4). 
Seventh, the same application brings to light certain practical advantages of the 
use of propensity scores which to our knowledge have not previously been noted 
(§ 4.3). 

1.2. A case study 

Our application is to Connors et al.'s (1996) highly influential, and contro- 
versial, study of the critical-care procedure known alternately as Swan-Ganz, 
pulmonary artery or right heart catheterization (RHC). RHC is a procedure to 
perform continuous measurements of blood pressure in the heart and large blood 
vessels of the lungs. Introduced in 1970, it became standard procedure without 
first being tested in clinical studies, as might be expected today, and empirical 
assessments that were subsequently conducted failed to uncover evidence that it 
improved medical outcomes (e.g., Gore et al., 1987; Zion et al., 1990). However, 
these studies were criticized for insufficient confounder controls. Using a large 
sample, good measures and extensive adjustments for background variables, 
Connors et al. (1996) echoed the disappointments of the earlier assessments and 
went further, finding RHC to worsen, rather than improve, both mortality and 
the duration of treatment in critical care. Each of these studies used nonrandom- 
ized data and is in some degree vulnerable to omitted variable bias. Although 
the results of subsequent randomized trials have been largely consistent with this 
picture (Rhodes et al., 2002; Sandhan ct al., 2003; Shah and Stevenson, 2004; 
Richard and et al, 2003; Harvey et al., 2005), the procedure remains a staple of 
critical care, and the surrounding debate continues. Our illustration examines 
how the omission of covariates from Connors et al.'s data might skew the point 
and interval estimates they support, specifically covariance-adjusted estimates 
of RHC's effects on length of hospital stay. 

1.3. The SUPPORT data 

Connors et al.'s (1996) data come from the Study to Understand Prognoses and 
Preferences for Outcomes and Risks of Treatments (SUPPORT). The study col- 
lected data on the decision-making and outcomes of seriously ill, hospitalized 
adult patients. Patients included in the study had to meet certain entry cri- 
teria and a predefined level of severity of illness. All 5735 SUPPORT patients 
who were admitted to an ICU in the first 24 hours of study were analyzed 
together, for reasons detailed in Connors et al. (1996). Data includes initial dis- 
ease categories upon admission, physiological measurements, and demographic 
information. Patients were coded as having had RHC if it was performed within 
the first 24 hours of the study. 
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Length-of-stay is the difference, in days, between a patient's recorded dates 
of entry and exit from the hospital. Some patients' stays ended in death, com- 
plicating interpretation of RHC's effects on length-of-stay: long stays are not 
in themselves desirable outcomes, but if RHC were to reduce mortality then 
its use might increase lengths of stay by keeping some patients in the hospital 
rather than the morgue. A full treatment of the issue is beyond the scope of 
this paper; we note, however, that Connors et al. found available data to speak 
against such a scenario, as does a supplemental analysis (for a subgroup with 
elevated mortality) that we present in Section 4.2. 

1. 4- First-pass regression results 

The outcome, length-of-stay, is right-skewed. We log-transform it before re- 
gression analysis. There are approximately 50 covariates for which regression 
adjustments might be considered; a backward stepwise selection procedure re- 
duces this number to 19, and estimates the RHC effect as 0.11: taking the log 
transformation into account, RHC seems to increase lengths of stay by about 
100 * (exp(.ll) — 1) = 12%. To reflect variability added by variable selection 
(Faraway, 1992), we ran a nonparamctric bootstrap to determine null quan- 
tiles for regression parameters' ^-statistics; in this case, however, bootstrap t- 
quantiles were similar to quantiles of the appropriate i-distribution. Either way, 
a 95% confidence interval for the RHC effect ranges from .06 to .16, encompass- 
ing increases of 6% up to 18%. 

2. Effect estimates accounting for omitted variables 

To understand how an omitted variable, W, could affect the coefficient on the 
treatment variable, Z, in a hypothetical regression of a given outcome on these 
and other variables, 

Y = aX. T + j3Z + SW + e, (1) 

it is well to begin by examining how included variables affect the treatment 
coefficient in the regression that was actually performed, 

Y = aX T + bZ + e. (2) 

This process lends context to an accompanying sensitivity analysis. 

2.1. Omitted variable bias 

Perhaps the most familiar way of comparing regressors is in terms of their effects 
on R 2 ; we begin in this mode. Of 33 regressors, the one that moves R 2 the most 
is "DNR status," an indicator of whether there was a do-not-resuscitate order 
on file for the subject on the day that he first qualified for inclusion in the 
study. Without it, R 2 would be reduced from .141 to .112; in contrast, removal 
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of the next-most predictive regressor reduces R 2 only to .131. On this basis 
one might expect DNR status to have a relatively large effect on inferences 
about the RHC effect, and in a sense it does: omitting it from the regression 
equation increases the RHC coefficient estimate from .112 to .143, 1.2 standard 
errors. For comparison, consider a regressors which contributes more modestly 
to the outcome regression: omission of "bleeding in the upper GI tract," for 
instance, reduces R 2 only by .001, the 28th smallest of 33 such reductions, and 
removing it moves the treatment coefficient only a few percent of an SE. Blurring 
the distinction between parameters and estimates just a little, we refer to this 
difference as bias, the bias that would have been incurred by omitting the DNR 
variable, ft appears that R 2 does track omitted-variable bias, at least to some 
degree. 

There is a simple relationship between a covariate W's contribution to ex- 
plaining the response, as measured by its impact on R 2 , and the bias incurred 
by its omission. According to Proposition 2.1, this bias is the product of the 
standard error on Z, the treatment variable, as calculated without regression 
adjustment for W; another term to be explained presently; and p y . w \ zyL , the par- 
tial correlation of W and the response given Z and remaining covariates, X . In 

turn > P 2 yu,\zx C£ l uals [( 1 - R low)-( 1 - R ldth.w)}/( 1 - R low)^ thc proportionate 
reduction in unexplained variance when W is added as a regressor (Christensen, 
1996, Ch. 6). That is to say, the square of the bias due to omitting a covari- 
ate is linear in the fraction by which that covariate would reduce unexplained 
variance. Returning to our example, DNR status consumes 3.3% of outcome 
variation that would otherwise be unexplained, which is 16 IQRs greater than 
the third quartile of Py. w i zx 's associated with the 33 available regressors; the 
Py-w\zx f° r u pper GI bleeding, in contrast, is .1%, just above the first quartile. 

The remaining factor in Proposition 2.1's expression for omitted variable 
bias expresses the degree to which the variable is confounded with treatment. 
This confounding turns out to figure more centrally in omitted variable bias, 
despite the fact that measurements of it arise less often in regression analysis 
than do p y . w \ zx or R 2 . A straightforward descriptive measurement falls out of 
the regression of Z on the remaining regressors from the outcome regression: 
their confoundedness with Z is reflected in their ^-statistics, the ratios of their 
coefficient estimates to their conventional standard errors. Note that these are 
here considered as descriptive, not inferential, statistics; were this auxiliary re- 
gression instead being used for statistical inference, an nonlinear model may be 
required for binary Z, in which case these nominal standard errors may not be 
appropriate. Denoting by tw the i-ratio attaching in this way to a regressor W, 
by b the RHC-coefficient in the absence of W and by (3 the same coefficient when 
W is included, the bias due to omission of a regressor decomposes as follows. 

Proposition 2.1. If R 2 . ZXL < 1 and tw is finite, then 

b- $= SE(b)t w p y .w\zx- ( 3 ) 
A proof is given in Section 3.2. 
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Besides being a strong predictor of the response, DNR status ranks highly 
in terms of its confounding with the treatment: if RHC is regressed on the 33 
covariates then its ^-statistic is —6.5, the fifth-largest in magnitude. Consistent 
with there being little bias associated with removing it, upper GI bleeding is 
as weakly related to the treatment variable as it was to the response; its t is 
only — .7, the fifth-smallest in magnitude. The strongest confounder among the 
adjustment variables is the Pa02/FI02 ratio, a reflection of oxygenation in the 
lungs, with t = —15.3. On the other hand, its relatively small contribution toward 
explaining variation in the outcome, Py. w u x = .0007, limits the effects of its 

removal on the treatment coefficient; here b—f3 is .011, or 40% of an SE. Although 
the effect on the treatment coefficient of excluding the Pa02/FI02 variable is 
tempered by its small corresponding Py, w i zx , (3) says that an actual omitted 
variable this confounded with treatment could contribute up to 15.3 standard 
errors' worth of bias. If this possibility is not merely notional, it follows that 
uncertainty assessments based in the usual way on SE(6) are quite incomplete. 

Proposition 2.1 suggests that the degree of an omitted variable's confounding 
with the treatment plays a potentially much larger role in determining the bias 
from its omission than does its conditional association with the response. As it 
can explain no more than 100% of otherwise unexplained variation, in general 
Py-w\zx — 1> ensurm g in turn that \b — /3|/SE(6) < \tw\- In contrast, nothing 
bounds the f-ratios tyy — covariates highly collinear with the treatment can 
introduce biases that are arbitrarily large, at least in principle. We refer to 
an unsigned statistic \t\y\ as the treatment confounding of W a covariate or 
potential covariate W. 

Table 1 reports tyy and p y . w \ ZXL values that actually occur in the SUPPORT 
data, placing observed covariates in the role of W one at a time. With these 
data, tw is often of a magnitude to make large biases possible (although this 
is somewhat tempered by relatively small p y . w \ Z x s )- The calculations permit 
the statistician and his audience to calibrate intuitions for tw against actual 
variables and their relations to the treatment. We refer to this as benchmarking 
of treatment-confounding. 

2.2. Omitted variables and the treatment's standard error 

The representation (3) of omitted-variable bias, as the product of SE(fo) with 
other factors, lends hope that a meaningful combinations of errors due to vari- 
able omissions and sampling variability might be bounded by a multiple of 
SE(6). A tempting solution is simply to add to the normal 97.5 percentile 1.96, 
or to whatever other multiplier of the SE has been found appropriate to de- 
scribe likely sampling error, a quantity large enough to bound plausible values 
of \twPy-w\zx\- This solution will often miss the mark, however, since omitting 
a covariate affects standard errors as well as coefficient estimates. This is illus- 
trated in Table 2, which shows effects on estimation of the RHC-coefficient of 
omitting one- by-one the regressors discussed in § 2.1. 

According to Table 2, adjusting the RHC coefficient for covariates includ- 
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Table 1 

Selected included covariates' relation to treatment and response variables, given remaining 

included variables 





L U111U L111LL111H, 


/0 tJ.L.t^l L^CtSL. Ill LllH_.A.]JlCtlllt_-Ll 




witil KriL, (Jtw / 


variation by adding W 




(rounded) 


( 1 DOn^ ^ 
' y-w\zyi ' 


income 


6.8 


0.3 


initial disease cat. 1 


48.1 


3.4 


initial disease cat. 2 


20.2 


0.8 


Comorbidities illness: 






renal disease 


2.1 


0.2 


Upper GI bleed 


0.7 


0.1 


Day 1 Measurements: 






APACHE score 


5.1 


0.1 


white blood cell ct 


0.5 


0.0 


heart rate 


2.5 


0.0 


temperature 


2.3 


0.1 


Pa02/FI02 


15.4 


0.1 


albumin 


2.3 


0.7 


hematocrit 


3.3 


0.9 


bilirubin 


2.2 


0.1 


sodium 


3.1 


0.1 


PaCo2 


6.8 


0.2 


DNR 


6.5 


3.3 


PH 


3.7 


0.3 


Admit Diagnosis Categories: 






neurology 


5.9 


0.2 


hematology 


3.6 


0.1 



ing DNR order gives a smaller standard error, .0260, than does excluding DNR 
order and adjusting for remaining covariates only (SE=.0264). In contrast, inclu- 
sion of PaO2/FI02 among the adjustment variables has the effect of increasing 
the standard error, from .0255 to .0260. Including or excluding upper GI tract 
bleeding leaves the standard error unchanged, suggesting that omitted variables 
associated with little bias should have little effect on the standard error. This 
turns out to be true, in a limited sense, as can be read from Proposition 2.2: the 
same statistics governing the bias due to omission of a covariate W also govern 
its effects on the standard error, in such a way that when both tw and p y . w \ zx 
are small then including or excluding W has little effect on the Z-cocfficient's 
standard error. 

Proposition 2.2. If R 2 y . zx < 1 and tw is finite, then 



SE0) = S E(b)^J±^-yJl-p y . w ^. (4) 

Here df=n — rank(X) — 1, the residual degrees of freedom after Y is regressed 
on X and Z . 

Proposition 2.2 will be strengthened in Proposition 4.3, which is proved in 
the Appendix. 
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Excluded 
Variable 


Estimate 


Standard Error 
df=5700 df = 50 


no exclusion 


0.112 


0.0260 


0.278 


DNR order 


0.143 


0.0264 


0.206 


GI bleed 


0.112 


0.0260 


0.274 


Pa02/FI02 


0.122 


0.0255 


0.115 



Table 2 

RHC- coefficient and its standard error after removing specified variables. 



Whereas \tw\ and and |/Cy.u,| zx both associate directly with the size of the 
bias from omitting W, they act in opposite directions on the standard error 
- reflecting the fact that increasing R 2 tends to reduce coefficients' standard 
errors, except when it's increased at the expense of introducing collinearity 
among regressors. The difference explains why omitting some variables increases 
standard errors whereas omitting others decreases them. 

It also reaffirms that omitted variables' effects on statistical inferences are 
only incompletely reflected in omitted variable bias. In the context of the present 
example, variable omissions have only modest effects on the standard error; but 
this is a consequence of the sample being quite large (n = 5735) relative to the 
rank of the covariate matrix (33). With df=5700, even a t\y of 16 inflates the 
standard error by at most 2%, according to (4). In moderate or small samples, 
standard errors will be more sensitive. The third column of Table 2 presents 
standard errors as they would be if the sample had been much smaller, such 
that df = 50, but with the same sample means and covarianccs. While the 
omission of a variable like GI bleed, which is only weakly related with already- 
included variables, still leaves the treatment's standard error much the same, 
variables that are either moderately or strongly confounded with the treatment 
now cause wide shifts in the magnitude of the standard error. Adjustment for the 
Pa02/FI02 measurement, for example, now increases the treatment's standard 
error by a whopping 59%. 

2.3. One sensitivity parameter or two? 

Should we be inclined to consider worst-case scenarios, Proposition 2.2 reaffirms 
Proposition 2.1's message that it's the omitted variable's treatment-confounding, 
not its potential to increase R 2 , that most demands our attention. The p y . w \ zx - 
contribution to SE(/3) is a factor bounded above by 1, whereas the ivK-factor can 
again be arbitrarily large — as could the ^-contribution to omitted-variable 
bias. The greater potential disturbance from confounding with the treatment 
than from large new contributions to the outcome regression seems to be a fea- 
ture multiple regression analysis shares with other statistical techniques; indeed, 
sensitivity analysis methods which set out to limit inference errors in terms of a 
single sensitivity parameter parametrize confounding with treatment, not pre- 
dictivity of the response, countenancing arbitrarily strong response-predictors 
in the role of the hypothetical omitted variable (e.g., Roscnbaum, 1988). 
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What would such an analysis suggest about the sensitivity of our RHC results 
to the omission of a confounder not unlike included variables? Pa02/FI02 was 
sufficiently confounded with treatment to indicate as many as 16 standard errors' 
worth of bias. Yet its inclusion in the regression adjustment moves the treatment 
effect by less than half a standard error. To have moved (3 so much it would 
have had to consume all of the variation in Y that was not explained without it, 
something no one would have expected it to do. To the contrary, inspection of 
other variables' contributions to i? 2 , as enabled by Table 1, suggests that .1 (for 
instance) would be a quite generous limit on \p y . w \ zx \. This in turn would restrict 
omitted variable bias due to a treatment-confounder as strong as Pa02/FI02 
to 1.6 SEs still a meaningful addition to the uncertainty estimate, if a less 
alarmist one. Rather than reducing the number of sensitivity parameters by 
permitting p y . w \ zx to vary within its a priori limits, it is more appealing to 
retain p y . w \ zx , restricting it within generous bounds of plausibility. 



3. Sensitivity Intervals 

Taken together, Propositions 2.1 and 2.2 enable a precise, closed-form descrip- 
tion of the union of interval estimates {/3 ± gSE(/3) : \tw\ < T, P y . w i zx < D }, for 
any nonnegative limits T, D on omitted variables' treatment-confounding and 
contributions to reducing unexplained variation — the collection of Z-slopes 
falling within the confidence interval after addition of a covariate W such that 
—T < tw < T and Py. w i zx < D. Such a union of intervals is itself an interval: 
following Roscnbaum (2002), we call it a sensitivity interval; and following Small 
(2007), we refer to the determining set of permissible values for (tw, p\. w \ zx ) as 
a sensitivity zone. The mapping of sensitivity zones to sensitivity intervals is 
given in Proposition 3.1. 

For p 2 £ [0, 1] and df w > define / by 



f(t,p 2 ,q,df w ) 



t\p\ + q^{l + ^-(\ + t 2 ))(\~ p 2 ), p 2 <g(t,q,df w ) 
^t 2 + q 2 -{l + ^-d + V)), p 2 > g(t, q, df w ); 

(5) 

where g(t,q,df w ) = sin 2 {arctan[tg- 1 (l + df w : (1 + £ 2 ))~ 1/2 ]}- 

Proposition 3.1. Let Y, X, Z and W be as in (1) and (2), with both regressions 
fit either by ordinary least squares or by weighted least squares with common 
weights. Let p y . w \ zx andtw be as defined in § 2.1; suppose Ry. zx < 1, \tw\ < oo- 
Fix q > 0. Then 



P±qSE0) Cb± f(\t w \,pl. wlzx ,df w ,q)SE(b) (6) 

where f is given by (5) and df w = df — 1 is the residual degrees of freedom in 
the regression ofY on X, Z, and W. 

Proposition 3.1 flows from the bias and variance expressions in Proposi- 
tions 2.1 and 2.2. Its proof appears in the Appendix. 
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3. 1 . Pegging the boundaries of the sensitivity zone 

Since our analysis began with covariate selection, there are a number of delib- 
erately omitted variables that can be used to peg at least one boundary of the 
sensitivity zone. As one might expect, the covariates eliminated by the stepwise 
procedure add little to those variables that were included in terms of prediction 
of the response, and it seems too optimistic to suppose of a genuinely unmea- 
sured confounder that its contribution to the outcome regression would be no 
greater than that of measured covariates that variable selection would put aside. 
On the other hand, confounding with Z plays little role in common stepwise 
procedures like the one we used, and the deliberately omitted variables can be 
used to guide intuitions about plausible values of tw- We selected six of these 
whose partial associations with RHC spanned the full range of such associations 
among stepwisc-climinatcd covariates and used their in/-valucs to delimit the 
first dimension of several sensitivity zones. Table 3 delineates the iiy-part of the 
sensitivity zone accordingly, choosing the bound T on treatment-confounding to 
coincide with the magnitude of confounding, conditional on stepwise-selected co- 
variates, between the treatment and each of our six covariates. Reference values 
for tw might also be obtained from external data sets. 

Table 3 

95% sensitivity intervals for the treatment coefficient, with the putative unobserved 
variable's treatment- confounding (\tyy\) hypothesized to be no greater than the 
treatment- confounding of 6 deliberately omitted variables. The decrease it would bring to the 
variance of response residuals is hypothesized to be no greater than either of 2 index values, 

1% and 10%, or is not restricted. 





Treatment 


% decrease 


in unexplained variation 




confounding 




(100p 2 , ) 




Variable 


benchmark 


1% 


10% 


Unrestricted 


Insurance class 


12.2 


most 


(0.03, 0.20) 


(-0.04, 0.26) 


(-0.21, 0.43) 


Respiratory eval. 


8.9 


some 


(0.04, 0.19) 


(-0.01, 0.23) 


(-0.12, 0.35) 


Mean blood press. 


8.6 


some 


(0.04, 0.19) 


(-0.01, 0.23) 


(-0.12, 0.34) 


Cardiovascular eval. 


8.5 


some 


(0.04, 0.19) 


(-0.01, 0.23) 


(-0.11, 0.34) 


Weight (kg) 


6.1 


some 


(0.04, 0.18) 


(0.01 , 0.21) 


(-0.05, 0.28) 


Immunosuppression 


0.4 


least 


(0.06, 0.16) 


(0.06 , 0.16) 


(0.06 , 0.16) 



Many analysts will have sharper intuition for potential covariates' effect on 
i? 2 , making it relatively easy to set plausible limits on Py. w i zx - Table 3 con- 
siders Py. w \ zx < -01 or .10, which may be useful as general starting points. In 
the present case study, for instance, when included covariates are removed and 
put in the role of W, p 2 y . w \ zyL — -01 corresponds approximately to the second 
most predictive of them, and the strongest single predictor (DNR order) gives 
Py-wlzx. = -03- P y . w \ zx < -1 appears to be quite a conservative bound: it's difficult 
even to find sets of included covariates that jointly contribute so much to the 
outcome regression as this. Only by simultaneously placing all of the covariates 
into the role of W, leaving the intercept alone in the role of X, does Py. w i zx 
reach .10. Restoring these variables into the regression, and barring scenarios of 
still-omitted variables for which Py. w i zx exceeds .1, Table 3 shows the result of 
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a positive treatment effect to be sensitive to omitted confounding on par with 
some of the strongest of the included confounders, but insensitive to confounding 
weaker than that. 

3.2. Basis for sensitivity formulas 

Propositions 2.1, 2.2, and 3.1 extend better-known descriptions of bias in re- 
gression coefficients' point estimates due to variable omission {e.g., Seber, 1977, 
p. 66) to interval estimates. They also have antecedents in earlier literature on 
numerical adjustment of multiple regression results for the addition or removal 
of a covariate (Cochran, 1938). Of the three, Proposition 2.1's proof is the most 
illuminating. It also conveys the flavor of the others, which appear in the ap- 
pendix. 

Consider X to be a matrix containing a column of l's (or columns from 
which a column of l's can be recovered as a linear combination) and let Y, 
Z, and W be column vectors of common extent, equal to the number of rows 
of X. An inner product is defined as (A, B) := J2 w i a il>i/ ^2w i: where lUj is a 
quadratic weight for the ith observation (in the case of unweighted least squares 
regression, Wi = 1). Write 1 for the column vector of Is. For vectors A, B, and C, 
let P](A\B, C) represent the projection of A into the subspace spanned by B and 
C. Variances and covariances are defined as follows: <r a b-c '■= {A — Pj(A\C), B — 
P](B\C)), u\. c = Oaa c^ab = &ab i, &a = °a i- Partial correlations are then given 
as: p ab := a ab /{a a a b ); p ab . c := o- ab . c / (o- a . c a b . c ). Denote the degrees of freedom 
available for estimating b as df = n — m — 1, where m = column. rank(JT) . The 
nominal standard error estimates for b and (3 (cf. (2) and (1)) arc then 

SE(S) =df~ 1/2 ^ and SE(/3) = (d/-l)- 1 / 2 ^^. (7) 

Proof of Proposition 2.1. To show b — (3 = SE(b)tw p y - w \zyn write 

Pj(W|Z,X) =: B*Z + C*X t . (8) 

Using (2), (1), and (8), the bias in the treatment coefficient can be written 
as b — (3 = B*S , a well-known result (Sober, 1977, p. 66). 

Write W ±x for W - Pj(VF|l,X), Z ±x for Z - Pj(Z|l,X), Y ±zx for Y - 
Pj(F|l,Z,X), and W-*-** forW-Pj(V|l,Z,X). Then Pj{W^ \Z^) = B*Z J ~*, 
and Pj(y- L -|W- L «) = SW^*. These formulas imply B* = a wz . x /a 2 z . x and 
o = o- yw . zx /a 2 w . zx = Pyw-xxVyzx/vw-zx, so that b - $ = B*5 can be written as 
the product of a y . zx /a z . x , (T wz . x /(a z . x a w . xx ) and p yw . zx . Introducing mutually 
cancelling factors of {df) ±x / 2 to the first and second of these and applying (7) 
turns this into the product of SE(6), (df) 1 ^ 2 a wz . x /(a z . x But 
tw is just the ratio of o wz . x j a 1 w . x to ov u , x /[(<i/') 1 / 2 c7. !i >-x], which simplifies to 
the second of these terms. The result follows. □ 
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4. Extensions 

As it is presented in Section 2, our method explores sensitivities of covariancc- 
adjusted estimates of a main effect to the omission of a single covariate. It 
may appear to be limited, then, to effect estimates made by linear covariate 
adjustment, without interaction terms or other allowances for heterogeneity of 
treatment effects, and to hidden bias scenarios involving omission of a single 
variable, rather than several. Such an appearance would be misleading. 

4-1- Several variables omitted at once 

Suppose now that W denotes not one but several omitted variables, or that it 
represents a single nominal variable with 3 or more levels, so that its encoding 
in terms of a design matrix would require 2 or more columns, and 2 or more 
degrees of freedom. Results previously presented still describe potential effects 
of W's omission, if tw is reinterpreted in a natural way. (p y . w \ zx retains its 
original interpretation, as the proportionate decline in unexplained variance 
from including W as a regressor.) 

When Z is regressed on X and a multivariate W, there is no one VK-coefHcient 
and corresponding i-statistic. The natural analogue of such a statistic is the 
ANOVA F-statistic comparing regression fits with and without W, F\y\ for 
univariate W, Fyy = t w , as is well-known. When rank(W) > 1, define the 
omitted variables' treatment confounding, again denoted tw, as the positive 
square root of [k(df)/(df + 1 — k)]Fw- Proposition 2.1 then gives the following. 

Corollary 4.1. Suppose R y . zx < 1, t 2 w is finite, and rank(W) = k > 1. Then 
{b - /3) 2 < V(b)[k(df)(df + 1 - ky^Fwp 2 or equivalent^ 



\b~P\ <SE(b)t W \ Py . wlzX \. 

Proof of Corollary 4-1- Without loss of generality, W is uncorrelated with Z 
and X: if not, replacing W with W — Pj(W\X 7 Z) leaves Z-coefficients and 
their standard errors unchanged. Define W = Pj(y ±x,z |M / ), where y ±x > z = 
Y - Pj(Y\X, Z). Again without loss of generality, W = (W, W 2 , ■ . . W k ), where 
W A. (W 2 ,...W k ). Writing 

Pj{Y\Z, X, W) =: a + PZ + jX T + 5{W + S 2 W 2 + ■■■ + S k W k , (9) 

it is immediate that 5 2 , ... ,5 k = , since W 2 , . . . , Wk are orthogonal toPj(F ±x ^|VK), 
and hence orthogonal to y- Lx ' 2 . Projecting (9) onto the span of Z, X, and then 
equating the .^-coefficient in what results with the Z-coefficient in (2) yields 

$ + S t Bi = b, (10) 

where Bf is defined by Pj(W\Z,X) = BfZ + C*X. In other words, b and (3 
are related just as they would have been had W been of rank 1, rather than k, 
consisting only of W. 
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We record some entailments of the definitions of p y . w \ zyi , tw and Fw in a 
lemma, proved in the appendix: 

Lemma 4.2. Suppose Ry. zx < 1, t w is finite, and rank(VK) = k. Then 

1 n 2 - n 2 9 t 2 < J? df Pw 

1 - Pyw\zx — PyS>\xx> a - l W K df+l-k rw 

The desired result now follows from (10), Proposition 2.1, and Lemma 4.2. □ 

When rank(W) > 1 we have the following variant of Proposition 2.2, proved 
in the appendix. 

Proposition 4.3. Suppose Ry. zx < I, t\f is finite, and rank(W) = k, k > 1. 
Then 

(i-pjUi-)- ( n ) 

Because the sensitivity intervals in Proposition 3.1 follow algebraically from 
the bias and standard error representations (3) and (4), they work for W of 
arbitrary rank. The proofs of Proposition 3.1 and the following are essentially 
the same. 

Proposition 4.4. fn the setting of Proposition 3.1 except with rank(14 7 ) > 1, 
(6) holds with df w =df- rank(W) and t w = {[k(df)/(df + 1 - k)]F w } 1/2 - 

Some of the hypothetical omissions discussed in Section 2 are of the type 
for which Proposition 4.4 is needed. The variable "Insurance class" appearing 
in Table 3, for example, is a nominal variable with 6 categories, consuming 
5 degrees of freedom when added as a regressor. Its treatment-confounding, 
tw, was calculated as the appropriately rescaled square-root of the f-statistic 
comparing the linear regression of Z on X and it to the regression of Z on X 
alone, about 2.24. 

4-2. Treatment effects differing by subgroup 

Recall from Section 2.1 that of the 33 X-variables selected as covariates for the 
regression of length-of-stay on RHC, DNR status most reduced R 2 . Patients 
with do-not-resuscitate orders suffered 25% greater mortality during the study 
than other patients, probably making it inevitable that their outcomes on this 
variable should systematically differ from patients without such orders. It is 
natural to suspect that effects of RHC might differ for them as well. In this 
case our linear models require interactions between RHC and DNR status - 
and perhaps other interactions as well, for that matter, but it suffices for us to 
restrict attention to the treatment interaction with a single binary moderating 
variable, as all issues pertaining to sensitivity analysis arise in this simplest case. 
Marcus (1997) explores related problems. 

Supplementing the regression of length-of-stay on covariates and RHC with 
an interaction between RHC and DNR status gave quite revealing results. The 



V0) = V(b) 



df + t 2 v 



df-k 
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additional right-hand-side variable, an indicator of RHC and DNR simultane- 
ously, bears a coefficient of —.43 and a ^-statistic of —5: it appears that the model 
devoting a single term to the treatment obscured significant effect heterogene- 
ity. Correspondingly the main RHC coefficient, interpretable as the effect for 
patients without DNR orders, is larger (+.15) than it was in earlier analyses 
without interactions (+.11); its standard error increases slightly, from .026 to 
.027. 

To subject these results to sensitivity analysis, we'll again use index values 
to limit the impact of the hypothesized omitted covariate on R 2 : Py. w \ zx < -10 
remains a generous limit, as simultaneously adding all 33 covariates to the re- 
gression of length of stay on RHC (now interacted with DNR status) decreases 
unexplained variation by only slightly more, about 11%. To set suitable lim- 
its on the omitted variable's treatment confounding, imagine for the moment 
that the inclusion of an interaction term had been handled somewhat differ- 
ently: rather than adding ZXi, the product of RHC and DNR indicators, add 
in X = ZX\ — ~P]{ZX\\X, Z). The coefficient of this variable lacks any straight- 
forward interpretation, but its addition to the right hand side of the equation 
has precisely the same effect on remaining coefficients as would the addition of 
ZXi itself. To benchmark treatment-confounding, we would then regress Z on 
X and X. By construction, however, X is orthogonal to Z and X, so that it 
itself earns a i-statistic of in this fitting, and its inclusion has no effect other 
than to remove a single residual degree of freedom. In other words, the tw~ 
benchmarks extracted by regressing Z on X alone, used for sensitivity analysis 
of the RHC effect in the absence of interactions, serve just as well here after of 
multiplication by the factor [(df + l)/(c?/)] 1 / 2 (which in this case is effectively 
1). The first 2 columns of Table 4 use benchmarks gathered in this way, finding 
the conclusion that RHC increases lengths of time in the hospital for patients 
without DNR orders to be a bit less sensitive to hidden bias than was the anal- 
ogous conclusion for the analysis assuming homogeneous treatment effects, in 
Table 3. 







Treatment confounding 


benchmarks found 


with estimated 








'feet on patic 


nts 


effect on patie 


nts 


ETT- 


weighted a^ 








/out DNR oi 


der 




^ith DNR ore 


er 




efiect 




Variable 


tw 


r y-w\ zx 


< 0.1 


t w 


P 2 1 
r y - w\zx 


< 0.1 


t W 


r y-w\ zx 


< 0.1 


Insurance Class 


12.2 


(-0.01, 


0.30) 


11.9 


(-0.74, 


0.17) 


12.2 


(-0.03, 


0.27) 


Respiratory eval. 


8.9 


(0.02 , 


0.27) 


8.1 


(-0.64, 


0.08) 


-4.1 


(0.03 , 


0.20) 


Mean blood press. 


-8.6 


(0.02 , 


0.27) 


-8.2 


(-0.64, 


0.08) 


-5.1 


(0.03 , 


0.21) 


Cardiovascular eval. 


-8.5 


(0.02 , 


0.27) 


-7.7 


(-0.63, 


0.07) 


-7.1 


(0.01 , 


0.23) 


Weight (kg) 


6.1 


(0.05 , 


0.25) 


6.4 


(-0.60, 


0.03) 


-5.3 


(0.03 , 


0.21) 


Immunosuppression 


0.4 


(0.09 , 


0.20) 


0.4 


(-0.44, 


-0.12) 


-5.9 


(0.02 , 


0.22) 



Table 4 



Sensitivity intervals for subgroup effects and for weighted average effects 

When we instructed it to include the RHC-DNR interaction among the ex- 
planatory variables, our software might equally well have added an indicator of 
RHC and the absence of DNR, Z(l — X{). In this case, the main effect would 
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be interprctable as the effect of RHC for patients with, rather than without. 
DNR orders. It follows that had that effect been the object of our interest, we 
could construct a sensitivity analysis for it in the same manner as just above, 
by persuading our regression program to expand the interaction differently. In 
actuality things are still simpler than that; we really only need to take ordinary 
care in interpreting the regression results, and somewhat modify the benchmark- 
ing equation. The effect for patients with DNR orders is the sum of the main 
RHC effect and the RHC-and-DNR effects, .148 + -.430 = -.28, with estimated 
variance equal to the sum of the two estimated variances and twice their covari- 
ance, .0065 = (.080) 2 . In parallel, to benchmark treatment confounding for this 
analysis, regress the sum of the main RHC indicator and the interaction term, 
the product of RHC and DNR indicators, on covariates. This gives somewhat 
different results than did the benchmarking for the RHC effect on patients with- 
out DNR orders, which omitted the interaction term from the left-hand side of 
its regression equation; compare the middle and left columns of Table 4. 

The same approach yields a sensitivity analysis for any target parameter rep- 
resentable as a linear combination of main effect and interaction terms. Take the 
effect of treatment on the treated, or ETT, parameter, considered from within 
a model permitting treatment effects to vary within specified subgroups. If the 
groups arc, for simplicity, patients with and without DNR orders, then since 
7% of patients receiving RHC had DNR orders the ETT can be represented as 
the main effect plus .07 times the DNR-RHC interaction effect, estimated as 
.148 + (.07) — .430 = .118, with corresponding standard error .026. For bench- 
marking, regress on covariates the RHC indicator plus .07 times the product of 
RHC and DNR indicators, with results as given in the rightmost two columns 
of Table 4. 

4-3. Propensity- adjusted estimates of treatment effects 

Regression adjusts between-group comparisons by attempting to remove ad- 
justment variables' contributions from the outcomes before comparing them. In 
contrast, adjustments based on propensity scores attempt to divide the sample 
into strata within which treatment and control subjects have similar distribu- 
tions of covariates. We estimated propensity scores using all 50 of the SUPPORT 
data's covariates in a logistic regression (Rosenbaum and Rubin, 1984), finding 
six equally-sized subclasses made treatment-control differences on the covariates 
jointly insignificant at level a = .10 (Hansen and Bowers, 2008). One can couple 
such a stratification with linear modeling to estimate treatment effects. In the 
simplest variant, responses are regressed on the treatment and fixed stratum 
effects. Fitting such a model to the SUPPORT data gives an RHC effect similar 
to what was estimated after ordinary covariance adjustment, as in Section 1.4, 
but with somewhat larger standard errors. 

The main assumption for this model is that so far as the outcome, length of 
stay, is concerned, the only systematic difference between the RHC and non- 
RHC patient in a propensity stratum is RHC itself. Relax this assumption in 
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Table 5 

Sensitivity intervals for the treatment effect after ordinary covariance and propensity- score 
adjustment, illustrating that propensity adjustment better limits sensitivity to the omission 
of adjustment variables. For covariance adjustment, tw is limited by the confoundedness 

with the treatment of 6 variables that had been eliminated by a preliminary 
variable-selection procedure, as in Table 3; for propensity adjustment, limits on treatment 
confounding are set by separately removing each of these and calculating their t\y 's after 
propensity adjustment for remaining variables. 
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Propensity adjusted 


regression 




\t m \ 






* < -1 


\t w 
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Insurance class 


12.2 


(0.03, 0.20) 


(-0.04 


, 0.26) 


8.6 


(0.02, 0.18) 


(-0.03, 0.23) 


Respiratory eval 


8.9 


(0.04 , 0.19) 


(-0.01 


, 0.23) 


3.1 


(0.04, 0.17) 


(0.02, 0.19) 


Mean blood press. 


8.6 


(0.04 , 0.19) 


(-0.01 


, 0.23) 


6.8 


(0.03, 0.18) 


(-0.01, 0.22) 


Cardiovascular eval 


8.5 


(0.04 , 0.19) 


(-0.01 


, 0.23) 


5.4 


(0.03, 0.17) 


(0, 0.20) 


Weight (kg) 


6.1 


(0.04, 0.18) 


(0.01 , 


0.21) 


5.1 


(0.03, 0.18) 


(0, 0.21) 


Immunosuppression 


0.4 


(0.06 , 0.16) 


(0.06 , 


0.16) 


0.5 


(0.04, 0.16) 


(0.04, 0.16) 



favor of another to the effect that so far as differences between outcomes and 
their projections onto an omitted variable are concerned, within strata RHC 
and non-RHC patients do not systematically differ. Were the omitted variable 
to become available, we could adjust by adding it to the explanatory side of the 
linear model. Without it, we can do sensitivity analysis. 

Benchmarking treatment-confounding levels takes a bit more effort than be- 
fore: rather than simply regressing Z on covariates and recording their t- or 
-F-statistics, we have to account in some way for the propensity stratification. 
We do this by removing the covariates, one at a time, from the propensity model, 
after each removal subclassifying the sample into sextiles, as before, but now 
using the modified propensity score; then regressing Z on the withheld covariate 
and on the propensity strata in order to associate a t- or ^-statistic with that 
covariate. Results of this process appear in Table 5. 

The results exhibit a striking pattern: adjustment based on propensity scores 
gives causal inferences that are far less sensitive to omitted variables than does 
regression-based covariate adjustment. With it, one can expect less residual 
confounding with the treatment than with covariate adjustment, as seen in 
smaller i^-values on the propensity-score side of the table. In a sense, propensity 
scores focus on confounding with the treatment, whereas covariate adjustment 
focuses on covariates and the response. Recall from Section 2 that while both 
matter to omitted variable bias, confounding with the treatment is both more 
difficult to pin down and potentially more pernicious. It stands to reason that 
while propensity adjustment may pay a slight penalty up front, in terms of 
somewhat larger standard errors than covariance adjustment, it offers a greater 
return downstream, in reduced sensitivity to hidden bias. 

5. Summary 

For effect estimates adjusted for covariates using ordinary least squares, im- 
pacts of covariate omission on point estimates and on standard errors have been 
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represented in terms of two statistics related the omitted variable to included 
ones, a measure of how adding the variable would affect R 2 and a measure of 
its association with the treatment variable given included variables. We refer to 
the latter as the omitted variable's treatment-confounding measurement. When 
generous limits on how the omitted variable would affect R 2 can be defended, 
they yield far less pessimistic assessments of sensitivity than would be possible 
without such a limit. Unlike the statistic pertaining to R 2 , plausible limits on the 
treatment confounding statistic are unlikely to emerge from intuition alone; on 
the other hand, it is straightforward and informative to determine study-specific 
benchmarks for it using available data. 

The changes to the treatment coefficient's point and error estimates that the 
addition of an omitted covariate would cause have been represented as multiples 
of its standard error. So these representations yield error appraisals accounting 
for certain hidden biases in familiar terms, as a multiple of the SE. The method 
adapts readily to scenarios of multivariate omission, heterogeneous treatment 
effect, and combinations of regression with propensity scores. 
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Appendix A: Proofs of theoretical results 

Proof of Lemma 4-2. Under the conditions of the lemma, (1) and (2) can be 
established: (1) In a regression of Z and X on Y, adding W has the same effect 
on the Z-cocfficient and model R 2 as adding W. Thus, (1) holds. 
(2) Furthermore, W e span(lU), so W explains no more variation in Z than does 

W. Thus, R\ Mx > R\ Mx , which implies j^f^ > m • From thc defini- 

tion of the ANOVA F-statistic: F w = [(a% x - <r,. xw )/(A;)]/[of. xw /(d/ + 1 — k)], 
which implies F W = [(df + 1 - k)/k][p 2 z , w{ J(l - 

In Section 4.1, tw is defined for multivariate W. Since k = rank(PF), in case 
of univariate W, t 2 w = F w . As rank(VF) = 1, = ^ = df ■ The 

result follows. 

□ 

Proof of Proposition 4-3. To relate SE(6) and SE(/3), begin by rewriting some 
of the variance terms: <r 2 y . zxw = Vax(Y Xxxw ) = Var(y ±zx -Pj(Y ±zx \W Xz *)) = 
o- 2 y ZX - °~ 2 yzxP 2 yw . zx , which implies a y . zxw = (T y . zic (l - p 2 yw . zx )~ 1/2 ■ 

Since p 2 . w ^ is the proportionate reduction in residual variance when W is 
added to the regression of Z on X, a 2 z . xw = cr 2 . x (l - p 2 . w][X ). 

Thus, SE(/3) = (df - k)^/ 2 (a y . zxw /a z . xw )[(l - p 2 yw . zx )/(l - p 2 zw . x )} 1/2 , or 

equivalently, SE(b)(df/df - fc) 1/2 [(l - p\ w .^)l(l - P 2 ZW . X )} 1/2 - 

Now, (1 - pI.^Y 1 = 1 + k(df + 1 - k^Fw. By the relationship of F w to 

tyy, it follows that (1 — p 2 . w \ x )^ 1 = (df + t 2 v )(df~ 1 ), which combines with the 
SE relationship to give (11). □ 

Proof of Proposition 3.1. There are two cases. First suppose p y . w \ zx > g(t w ,q, c). 
For a fixed q and t w , by Proposition 2.1 and Proposition 2.2, the desired confi- 
dence interval $±qSE(0) can be represented as $— gSE(/3) = 6+Z(arcsinp J ,., i ,|j. x ) 
and (3 + gSE(/3) = b + u(axcsmp y . w i zx ), where arcsinp y . TO | zx G (— 7r/2,7r/2) and 
1(9) := asinf? — dcos9 and u(9) := asm9 + dcos9 for a := — twSE(b), and 
d := ((ty/ + df)/(df — 1)SE(6)) 1 ' 2 . By calculus, u(-) can be seen to have its max- 
imum over (— 7r/2,7r/2) at arctan (a/d) E (— tt/2, tt/2). By algebra and trigono- 
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metric identities, that maximum is (a 2 +d 2 ) 1/2 . Similarly, over (— 7r/2, 7t/2), l(-) 
takes its minimum value of —(a 2 + d 2 ) 1 / 2 at arctan — (a/d). 
Now suppose pj^^ > g(t w , q, c). For some p such that 

< p < sin arctan(|i w |[(i^ + d f)/( d f - 1)] _1/2 ), let /0 y .„,| 2X £ Then, 
argmaxp y iu|ax ^(arcsinpy.^izx) = p if <yt" < 0, otherwise it is — p if > 0. 
Similarly, argmin PH , x Z(arcsin/5 y . !U | 2X ) = — p if t\\* < 0, otherwise it is p if 
> 0. 

To verify the relationship for u. note w(arcsin(-)) takes only one maximum as 
its argument ranges over [—1,1] at sin arctan \tw\ [(t\v + df)/(df — l)] -1 / 2 > p. 
Thus, u(arcsinp. y . u ,| zx ) must take its maximum as p y . w \ zx ranges over [— p, p] at 
one of the endpoints of that interval. The analogue for I is established by similar 
reasoning. □ 
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